A Model-Based Framework to Identify Optimal Administration Protocols for Immunotherapies in Castration-Resistance Prostate Cancer

Simple Summary Although ipilimumab has been approved for the treatment of many types of cancer, most prostate cancer patients seem to not respond well to the therapy. Here, we mathematically investigate the clinical relevance of ipilimumab, as mono-therapy and in combination with the dendritic cell vaccine sipuleucel-T, for the treatment of castration-resistant prostate cancer patients. The employed optimization problem, which incorporates a function of toxicity depending on the drug concentration, establishes precise protocols of administration of ipilimumab to control or eradicate prostate tumor, and defines how changing of key parameters affects the outcome. Overall, this mathematical study can help in optimizing the clinical use of ipilimumab for the effective treatment of castration-resistant prostate cancer patients. Abstract Prostate cancer (PCa) is one of the most frequent cancer in male population. Androgen deprivation therapy is the first-line strategy for the metastatic stage of the disease, but, inevitably, PCa develops resistance to castration (CRPC), becoming incurable. In recent years, clinical trials are testing the efficacy of anti-CTLA4 on CRPC. However, this tumor seems to be resistant to immunotherapies that are very effective in other types of cancers, and, so far, only the dendritic cell vaccine sipuleucel-T has been approved. In this work, we employ a mathematical model of CRPC to determine the optimal administration protocol of ipilimumab, a particular anti-CTLA4, as single treatment or in combination with the sipuleucel-T, by considering both the effect on tumor population and the drug toxicity. To this end, we first introduce a dose-depending function of toxicity, estimated from experimental data, then we define two different optimization problems. We show the results obtained by imposing different constraints, and how these change by varying drug efficacy. Our results suggest administration of high-doses for a brief period, which is predicted to be more efficient than solutions with prolonged low-doses. The model also highlights a synergy between ipilimumab and sipuleucel-T, which leads to a better tumor control with lower doses of ipilimumab. Finally, tumor eradication is also conceivable, but it depends on patient-specific parameters.


Introduction
Prostate cancer (PCa) is one of the most common tumor in male population [1]. Although organ-confined tumors are successfully controlled by surgery or radiotherapy [2], Side effects due to drug administration are often classified as moderate, i.e., grade 1 and 2, and severe, with grade ≥3. The latter are generally more intense, and they can have a longer resolution time. The AEs observed in a clinical trial are always collected and analyzed in order to estimate the drug toxicity. Sometimes, promising clinical trials have been interrupted due to unexpected side effects, and therefore, this information cannot be ignored. In the literature, there are several examples of PK/PD mathematical models defining the drug toxicity also from a molecular point of view. For example, the work by El-Masri et al. [16] describes the toxicity by considering the effect of a mixture of drugs on the human body. Similar methods have been employed also in the context of prostate cancer [17]. Even if this approach have been previously used to define optimal control problems [18], for our scope, we do not need such a detailed description, which would increase the model complexity both from mathematical and computational point of view. Another possible approach for taking into account drug toxicity in optimal control problems is based on connecting the toxicity with the drug concentration, by defining a maximum acceptable drug exposure and a critical concentration threshold, over which the drug becomes not tolerated [19,20]. However, the ipilimumab has not yet been approved for PCa and there are not enough information about maximal tolerated doses and drug exposure. Therefore, we introduce a toxicity function depending on ipilimumab concentration.
Up to our knowledge, in literature there are few mathematical models describing drug toxicity function, and the most detailed has been introduced by Hadjiandreou et al. [21]. The authors described the toxicity as a linear function of the drug concentration multiplied by a coefficient, representing the side effect magnitude. This coefficient has been defined by considering the side effects registered into the patient population, weighted by their severity. However, often adverse events can be considered negligible for low drug doses, while they become more frequent and strong in case of high doses, showing a non-linear relation with drug concentration. Therefore, we use a non-linear toxicity to fit toxicity coefficients, computed from literature experimental observations. For a given clinical protocol, the corresponding toxicity coefficient is defined by the percentage of patients showing moderate and severe adverse events, weighted by a function of the average resolution times of the side effects. The occurrence of side effects can produce a stop of the treatment affecting its efficacy. The need of side effects recovery leads to either a delay or a definitive stop of the treatment delivery. Thus the evaluation of the side effects recovery time is essential in evaluating the drug toxicity. Indeed, in patient populations with equivalent adverse events and frequencies, the life quality of those patients having lower average resolution times considerably improves, and therefore we can take the resolution time as an indication of a a lower drug toxicity.
The paper is organized as follows. In Section 2, we describe the experimental data used to estimate the toxicity function, we present the mathematical model and some technical aspects related to the optimization. We show the model results in Section 3, including the estimates of the toxicity function and the proposed optimal administration protocols, obtained by fixing different constraints. In Appendix B, we also analyze the synergy between the ipilimumab and sipuleucel-T, by means of a modified Bliss combination index, in order to further investigate how the two immunotherapies work in combination. In the Discussion we summarize the main contributions, linking the mathematical results with biological observations.

Experimental Data
In order to evaluate the toxicity of the drug ipilimumab, we consider experimental data from literature. We select the only two clinical trials on PCa patients involving a single infusion of ipilimumab [22] and a dose-escalating administration [23], in order to compare different administration protocols. The study by Slovin et  with external-beam radiotherapy, but we select only the data regarding the ipilimumab administered as mono-therapy. The study by Small et al. considers patients treated with a single dose of ipilimumab. The authors collected the frequencies of the single side effects and observed AEs of grade 3 in 2 over 14 patients. This study does not report the percentages of AEs of grade 1 and 2, but, comparing all these data, we are able to derive a range of patients presenting moderate AEs. All these information are summarized in Table 1. Table 1. Percentages of AEs registered in prostate cancer patients subject to ipilimumab therapy [22,23]. The two columns Grade 1-2 AEs and Grade 3-4 AEs represent the magnitude of the side effects, which is grade 1 and 2, and grade 3 and 4, respectively. The data are expressed in terms of number of patients and corresponding percentages (p M and p S in the Equation (2)). In Small et al., severe AEs were only of grade 3, and their number is clearly reported, while we infer a possible range of patients presenting AEs of grade 1 and 2.

Grade 1-2 AEs
Grade 3-4 AEs n (%) n (%) Slovin  In our analysis, the resolution times of the symptoms is an essential information to evaluate the drug toxicity. However, no one of the selected studies refers the average resolution times of the side effects. To extrapolate this information, we consider the studies on melanoma, which is one of the cancer for which ipilimumuab has been approved. In these studies, the average resolution times have been registered only in case of immune-related AEs, and they change depending on the sample. In melanoma patients immune-related AEs are the most common side effects of ipilimumab treatments [24], and this is confirmed for PCa patients [23]. Therefore, we selected the studies from Weber et al. [25] and Hodi et al. [26], both administering 3 mg/kg of ipilimumab every 3 weeks for 4 cycles, and we set the resolution times of moderate and severe AEs as an average between the ones registered in these studies, obtaining t M = 38.15 days and t S = 45.15 days, respectively.

Mathematical Model
We employ a mathematical model of tumor growth, interacting with the immune system, developed in our previous work [15]. The original model is composed by 8 ordinary differential equations (ODEs), describing the two forms of PCa, (castration sensitive and resistant, previously named androgen dependent and independent) and some key variables of the immune system, such as dendritic cells, Cytotoxic T cells (CTL), Regulatory T cells (Treg), the interleukin-2 (IL-2), as well as the evolution of the Prostate Specific Antigen (PSA) as linear combination of the PCa cells. The model also includes three treatments: the androgen deprivation therapy, the dendritic cell vaccine sipuleucel-T and the anti-CTLA4 ipilimumab. Since we aim at evaluating the ipilimumab efficacy on castration-resistant prostate cancer, we set the androgen deprivation therapy as the mainstay treatment, such that the androgen dependent PCa form can be ignored (CSPC ≈ 0). Moreover, by imposing the quasi-steady state approximation on the variable IL-2, we can further reduce the model, obtaining the following ODE system: where T, C, R, D and I p represent, respectively, CRPC cells, CTLs, Tregs, dendritic cells and ipilimumab concentration, while I L (T, C) = i 0 + e C T s+T describes the dynamics of IL-2, obtained by the quasi-steady state approximation. Table A4 in the Appendix D summarizes the parameter estimates. Figure 1 shows a schematic representation of the modelvariable interactions. The CRPC cells are assumed to grow according to a logistic term, while they can be killed by CTLs. The effect of the drug ipilimumab has been described by an additional term which increases the CTL tumor-killing capacity (the last term into the first equation). The CTLs are activated by the dendritic cells, and can also increase in number as effect of the clonal expansions due to IL-2. Treg cells repress the immune reaction by reducing the number of activated T cells. These immune cells are activated either by dendritic cells or IL-2. The dendritic cell vaccine has been described by an infusion of dendritic cells.
The starting point of all the model simulations has been set according to the initial condition of the patients in the Slovin et al. study. The authors considered a sample of patients previously treated with androgen deprivation therapy, that have developed to the castration-resistant state. To reproduce the initial condition of this sample, we simulate the complete system in [14] by imposing androgen deprivation therapy, and we find the timet such that PSA(t) corresponds to the average baseline value recorded in the patient population. Then, we choose as initial conditions for (1) the values of the variables of that system at timet, namely:

Computation of Toxicity Coefficients
In order to quantify the toxicity of a given protocol, on the basis of clinical observations, we introduce the coefficients associated to moderate and severe AEs as: where C M and C S represent the coefficient of toxicity of moderate and severe AEs, t M = 38.15 and t S = 45.15 are the average resolution times of the symptoms expressed in days computed as described in Section 2.1, g M and g S are the severity grades and p M and p S are the percentages of patients presenting the side effects. Table 1 provides the values of p M and p S with respect to the severity grades (g M = 1.5 and g S = 3.5), in different administration protocols. In this representation, adverse effects with resolution time of 1 day do not contribute to the toxicity coefficient. Higher resolution times increase the value of the term 1 − 1 t k , so that the corresponding C k grows. Finally, we define the toxicity coefficient of a given administration protocol j as: whereT ox represents the maximal toxicity value, i.e., the coefficient computed in case of all patients having severe AEs, with very high average resolution time such that 1 t k ≈ 0. Therefore, the toxicity coefficients are values between 0 and 1.
By considering the experimental data described in Section 2.1, and by assuming that the times for resolution of AEs are comparable between prostate cancer and melanoma patients, we can compute the toxicity coefficients corresponding to these administration protocols. The results are reported in Table 2. Table 2. Toxicity coefficients computed by considering the data reported in Table 1. C M and C S represent the coefficients of toxicity of moderate and severe adverse events employed to evaluate the toxicity coefficient (C Tox ) for a given protocol by the Equation (3). The coefficient of the single ipilimumab dose is represented by its average value ± standard deviation, since the corresponding reported data allow us to only estimate a possible range of values.

Toxicity Function
The aim of this work is to find an optimal protocol for ipilimumab administration by balancing the efficacy of a therapy in reducing tumor and the drug toxicity. To this end, we need to define a toxicity function for any feasible protocol of ipilimumab administration. We assume that the toxicity is a function depending on the drug concentration all over the time of the treatment, namely: The shape of the function depends on the parameters a, b and n that are determined by fitting the function (4) to the experimental data, as described in the Appendix A. In any case, the toxicity is 0 if no ipilimumab is administered, increases as the doses are larger and has a bounded value for any possible protocol, consistently with the definition of toxicity coefficient.

Optimization Problems
We define optimal administration protocols for ipilimumab those that minimize a weighted combination of average tumor size and overall toxicity over an appropriate time interval, that we choose as 5 years. We computed the optimal protocols, both for the case of administration of ipilimumab in mono-therapy and in combination with sipuleucel-T.
The optimal controls have been computed under several constraints, which have been suggested by clinical practice, and have also the effect of reducing the computational burden. Precisely, we assumed that the therapy period cannot be longer than 6 months and that there is a minimal interval between treatments, either 1 week or 3 weeks (the typical interval in clinical trials).
The minimal (0.3 mg/kg) and maximal (10 mg/kg) doses have been fixed according to the ipilimumab doses tested by clinical trials [23,27,28]. Moreover, we assumed that the ipilimumab concentration has to remain below a threshold I MAX p , for all the treatment period. The threshold I MAX p has been established on the basis of the maximal-dose clinical trial, i.e. one infusion of 10 mg/kg every 3 weeks for four cycles. Simulating the system for a patient of 88.5 kg (average patient weight according to Coletti et al. approximation [15]), we found the maximal concentration reached and set it as the threshold I MAX p = 1266.55 mg. Mathematically, finding the optimal protocol consists in solving the minimum problem defined by the following equation: where t f = 5 years, and the weights w 1 = 1 and w 2 = 500 have been empirically fixed in order to make comparable the two terms of the objective function O 1 =  Table 1). This problem consists in a combinatorial optimization, which is well suited for genetic algorithms. Therefore, the solutions have been obtained by the Genetic Algorithm of Matlab [29,30], implementing a MustiStart approach, in order to explore all the parameter space. For each problem, we find an optimal administration protocol for the ipilimumab, administered as mono-therapy or in combination with sipuleucel-T, administered by following FDA suggestions, i.e., three infusions of 50 million of cells, one every 2 weeks. The minimum problems are repeated by imposing different minimum time intervals between two infusions, i.e., ∆ t = 1 week and ∆ t = 3 weeks, as discussed above.
We also consider a Second Optimization Problem (OP2), in which we add the constraint that T(T f ) = 0, i.e., that the tumor is eradicated at the end of the treatment. Note that the solutions of (1) mathematically are never exactly equal to 0. However since T(t) in (1) is measured in billions of cells, we set it equal to 0 att if T(t) < 10 −9 .

Toxicity Function Estimates
The toxicity function defined in Equation (4) needs to be calibrated. Figure 2 shows the fitting between function outputs and toxicity coefficients, which has been obtained by following the procedure presented in Section 2.4. The chart shows with red filled points the values of the ipilimumab toxicity coefficients listed in Table 2, with the relative error bar in case of the single-infusion. The blue empty points represent the values of toxicity predicted by the toxicity function. The data points represent different administration procedures: the point with error bar is the toxicity coefficient in case of a single infusion of ipilimumab, while the others refer to intermittent protocols with different doses. Therefore, in order to easily compare the data fitting, we put on x-axis the maximal concentration of the drug ipilimumab. Comparison between toxicity coefficients (red filled points) and the outputs of the toxicity function (blue empty points) evaluated reproducing the same clinical protocols described in Table 3.
The point related to the single infusion is represented with the relative error-bar.
Overall, our function provides a good qualitative fit to the data. In particular, the toxicity function reproduces quite well the coefficient values in case of intermittent protocols, in which the error is between 0.04 and 0.08. The worst fit has been obtained in case of single-infusion (the point including error-bar), with an error of 0.1.
The estimates of the parameters in Equation (4) are reported in Table 3.

Optimal Administration Protocols
As described in Section 2.5, we defined two optimization problems. The first one (OP1) finds the less-toxic anti-CTLA4 schedule able to reduce tumor growth. Figures 3 and 4 show the results obtained by considering, respectively, ∆ t = 1 week or ∆ t = 3 weeks, as the minimum time interval between two infusions. The charts compare the outcomes in case of mono-therapy with ipilimumab (AC) and combination therapy with ipilimumab and sipuleucel-T (AC + V). These numerical results are summarized in Table 4. finds the less-toxic anti-CTLA4 schedule able to reduce tumor growth. Figures 3 and 4 show 232 the results obtained by considering, respectively, ∆ t = 1 week or ∆ t = 3 weeks, as the minimum 233 time interval between two infusions. The charts compare the outcomes in case of mono-therapy 234 with ipilimumab (AC) and combination therapy with ipilimumab and sipuleucel-T (AC+V).

235
These numerical results are summarized in Table 4. In all cases, the optimal administration  remains under the threshold. This implies that, in case of mono-therapy, when ∆ t = 1 week, the 240 optimal administration times are at weeks t 1 = 0, t 2 = 1, t 3 = 3 and t 4 = 4 ( Figure 3), while, 241 when ∆ t = 3 weeks, the optimal administration times are at weeks t 1 = 0, t 2 = 3 and t 3 = 6 242 ( Figure 4). Changing the minimum time between drug administrations, the model suggests 243 two approaches of drug administration: one with more frequent low doses infusions ( Figure   244 3), the other with fewer infusions of high doses (Figure 4). Compared to the mono-therapy, the 245 combination therapy AC+V is predicted to reduce more rapidly the tumor mass, involving a lower 246 dosage of ipilimumab, and thus causing a lower toxicity.

247
The other optimization problem OP2 finds the optimal protocol able to eradicate the tumor.    Table 4. Numerical results for the optimization problem OP1. Anti-CTLA4 has been administered as mono-therapy (AC) and in combination with dendritic cell vaccine (AC + V). Both the immunotherapies have been tested by taking into account the general constraints defined in Section 2.5 and an additional constraint concerning the minimum time interval between two infusions (∆ t ) fixed as 1 or 3 weeks. Column three and four list the infusion times and doses predicted by the model as the optimal schedules. The values O 1 and O 2 represent the two components of the optimization function, .
In all cases, the optimal administration protocols are predicted to control prostate cancer over the therapy period. The optimal schedules strictly depend on the constraint of the minimum time interval ∆ t . Indeed, the model suggests to administer high doses of ipilimumab as soon as possible, provided that the drug concentration remains under the threshold. This implies that, in case of mono-therapy, when ∆ t = 1 week, the optimal administration times are at weeks t 1 = 0, t 2 = 1, t 3 = 3 and t 4 = 4 ( Figure 3), while, when ∆ t = 3 weeks, the optimal administration times are at weeks t 1 = 0, t 2 = 3 and t 3 = 6 ( Figure 4). Changing the minimum time between drug administrations, the model suggests two approaches of drug administration: one with more frequent low doses infusions ( Figure  3), the other with fewer infusions of high doses ( Figure 4). Compared to the mono-therapy, the combination therapy AC + V is predicted to reduce more rapidly the tumor mass, involving a lower dosage of ipilimumab, and thus causing a lower toxicity.
The other optimization problem OP2 finds the optimal protocol able to eradicate the tumor. Figures 5 and 6 show the model suggested administration protocols in case of two different constraints about the minimum time interval between two infusions of either ∆ t = 1 week or ∆ t = 3 weeks, respectively. The charts compare the mono-therapy with anti-CTLA4 (AC) with the effect of the combination therapy with anti-CTLA4 and vaccine (AC + V). The results are presented in Table 5.
As observed for the other optimization problem, the optimal protocol depends on the fixed minimum time interval between two infusions (∆ t ), and the model suggests to administer anti-CTLA4 as soon as possible. In this case, some doses are slightly increased reltively to what shown in Table 4, and an additional infusion is necessary in order to eradicate the tumor. The results with the combination therapy confirm the possibility of tumor eradication by using a reduced amount of ipilimumab.    Table 5. Numerical results for the optimization problem OP2. Anti-CTLA4 has been administered as mono-therapy (AC) and in combination with dendritic cell vaccine (AC + V). Both the immunotherapies have been tested by taking into account the general constraints defined in Section 2.5 and an additional constraint concerning the minimum time interval between two infusions (∆ t ) fixed as 1 or 3 weeks. Column three and four list the infusion times and doses predicted by the model as the optimal schedules. The values O 1 and O 2 represent the first two components of the optimization function, namely O 1 = t f t 1 Tum(t)dt and O 2 = Tox(t f , I p (t)).

Sensitivity of Results to Parameter Estimates
The computed results are influenced by the model parameter estimates. In particular, the value of k AC , representing the maximal killing rate of tumor by CTL due to the drug ipilimumab, strongly affects the efficacy of the proposed administration protocols. The value of this parameter has been estimated in our previous work [15], by fitting experimental data referring to those prostate cancer patients having the greatest benefits from the ipilimumab therapy. Therefore, the value estimated for k AC is higher than what adequate for the average of the prostate cancer patient population, and all the model results represent the optimal case of well-responding patients. In order to assess the sensitivity of the results to the value of k AC , we analyze the model dynamics for lower values of this parameter. Fixing the optimal administration protocol for the mono-therapy showed in Figure 6, we perform a sensitivity analysis on k AC . Figure 7 shows how the model outcomes change as the parameter value is reduced.
These results highlight that by halving k AC , even if the administration protocol is able to initially reduce the tumor, the prostate cancer starts increasing again 2 years after the therapy start (orange line). If the parameter value is further reduced, the same therapy is not able to control the tumor growth (yellow line).
Given the strong influence of the parameter k AC on the results, we investigate if tumor eradication could be possible if its value is half the reference value (Table A4). Therefore, we repeat the optimization procedures presented in Section 2.5 with the reduced value of k AC . Figure 8 shows the model suggested optimal protocols in case of mono-therapy with ∆ t = 1 week. The figure compares the optimal administration protocols by considering the estimate value k AC with the one obtained with k AC /2.
We repeated the optimization procedure also by setting ∆ t = 3 weeks (see Appendix C, Figure A1) and for the combination therapy with ipilimumab and sipuleucel-T, with both ∆ t = 1 week ( Figure A2) and ∆ t = 3 weeks ( Figure A3). population, and all the model results represent the optimal case of well-responding patients. In 267 order to assess the sensitivity of the results to the value of k AC , we analyze the model dynamics for 268 lower values of this parameter. Fixing the optimal administration protocol for the mono-therapy 269 showed in Figure 6, we perform a sensitivity analysis on k AC . Figure 7 shows how the model 270 outcomes change as the parameter value is reduced. These results highlight that by halving k AC ,  (Table A4). Therefore, we 276 repeat the optimization procedures presented in section 2.5 with the reduced value of k AC . Figure   277 8 shows the model suggested optimal protocols in case of mono-therapy with ∆ t = 1 week. The  Figure A1) and for the combination therapy with ipilimumab and sipuleucel-T, with both 281 ∆ t = 1 week ( Figure A2) and ∆ t = 3 weeks ( Figure A3). 282 The model predicts that the tumor eradication is possible even if the value of the parameter 283 Figure 8. Comparison between administration protocols obtained with different values of the parameter k AC . The dynamics refer to the ipilimumab mono-therapy as result of the optimization problem OP2, with ∆ t = 1 week. Blue lines represent the model outcomes with the estimate value of k AC (the same presented in Figure 5), while orange lines represent the one with k AC 2 . The ipilimumab threshold has been highlighted with a black dotted line. The model predicts that the tumor eradication is possible even if the value of the parameter k AC is half the reference. However, the amount of administered anti-CTLA4 drug is substantially increased, as well as the length of the therapy. By comparing Figures 8 and A1, we observe that by fixing ∆ t = 1 week the number of ipilimumab infusions needed to eradicate the tumor is higher than the one suggested with ∆ t = 3 weeks. In the first case, the model suggests to administer as much drug as possible, one dose per week for the first two months. After this period, the administrations could be a bit relaxed. On the other hand, by imposing the constraint of ∆ t = 3 weeks, the predicted optimal administration protocol is constituted by two standard protocols (one infusion per week for 4 cycles) repeated with a one month break.
Similarly to the previous cases, when the dendritic cell vaccine is combined with the anti-CTLA4, the model suggests to administer a lower amount of ipilimumab compared to the mono-therapy optimal protocols (Figures A2 and A3). In particular, when ∆ t = 1 week, the administered doses of anti-CTLA4 needed to eradicate the tumor passes from 12 for the mono-therapy to 9 for the combination therapy. The numerical results for the optimal administration protocols in case of reduced k AC are listed in Table 6. By fixing the parameter k AC at 1/4 of the estimated value reported in Table A4, the model predicts that tumor eradication is no longer feasible with 6-month therapies. Table 6. Numerical results for the optimization problem OP2, with a halved value of k AC . Anti-CTLA4 has been administered as mono-therapy (AC) and in combination with dendritic cell vaccine (AC + V). Both the immunotherapies have been tested by taking into account the general constraints defined in Section 2.5 and an additional constraint concerning the minimum time interval between two infusions (∆ t ) fixed as 1 or 3 weeks. Column three and four list the infusion times and doses predicted by the model as the optimal schedules. The values O 1 and O 2 represent the first two components of the optimization function, namely O 1 = t f t 1 Tum(t)dt and O 2 = Tox(t f , I p (t)). To understand how the parameter uncertainty affects model results, we investigate how the proposed administration protocols vary by perturbing the parameter of the toxicity function. To this end, we implement the optimal control problem by increasing or reducing the parameter a of the Equation (4). Our results show that the model predictions are not influenced by the parameter uncertainty, since the optimal control outcomes remain identical after moderate perturbations of a parameter. The predicted optimal administration protocol changes only if the value of a is multiplied at least by a factor 5 or divided at least by a factor 2. Even in those cases, the model suggestions remain in line with the ones obtained with the original parameters estimates, proposing always high doses over a short time period, instead of a low-dose and prolonged therapy. Note that changing the weights in the optimization objectives (5) is equivalent to multiplying or dividing the toxicity function. Thus these results can be viewed also as a sensitivity analysis on the weights.

Discussion
In this work, we employ a mathematical model of human prostate cancer in order to determine optimal ipilimumab administration protocols able to reduce/eradicate tumor by balancing the treatment ability in reducing tumor and the drug toxicity. Ipilimumab is an anti-CTLA4 approved for the treatment of several tumors, and tested in metastatic castration-resistant prostate cancer [10]. We consider patients previously treated with androgen deprivation therapy, who develop the castration-resistant prostate cancer form, and we investigate the efficacy of immunotherapies on those tumor cells. We examine the effect of ipilimumab administered as mono-therapy and in combination with the dendritic cell vaccine sipuleucel-T, administered following the FDA recommendations. Our results highlight that the administration of ipilimumab is potentially able to control or eradicate the tumor. In particular, the optimal administration protocols seem to be feasible, and the corresponding toxicity profile is comparable with those observed in the clinical trials (Tables 2, 4 and 5), suggesting that the proposed therapy could be well-tolerated. Moreover, the result obtained by fixing ∆ t = 3 weeks for the optimal administration of ipilimumab as mono-therapy can be compared to the high-dose tested protocol [23,28].
The combination with the vaccine is predicted to improve the efficacy of the monotherapy, since the corresponding optimal protocol provides a faster tumor reduction, while administering a lower amount of anti-CTLA4, with a consequent reduction in ipilimumabrelated toxicity. This effect is probably due to the increase in CTLs proliferation induced by the dendritic cell vaccine [31], which, coupled with a drug aiming at increasing the CTL tumor-killing activity, causes a stronger tumor suppression. However, clinical studies highlight that only few patients show an evident positive effect as a consequence of ipilimumab therapies. In the study by Small et al., only 2 out of 14 patients treated with a single dose of ipilimumab had a significant PSA reduction (>50%), while Slovin et al. [23] observed a complete and a partial response in 2 out of 16 patients, treated with intermittent protocols of high-doses. These different outcomes can be explained by patient heterogeneity, since the success of a therapy depends on patient-specific parameters. One of these is the parameter k AC , which measures how much the ipilimumab can improve the tumor-killing activity of CTLs. As shown in Figure 7, a perturbation of this parameter can qualitatively change the model outcomes, and the optimal protocol becomes unable to even control the tumor growth.
To further investigate if the tumor eradication could be possible even for those patients who respond less well to the standard therapies, we find the optimal protocol for lower values of k AC . By analyzing the optimal administration protocols we obtained under different model constraints on ∆ t , i.e., the minimum time interval between two infusions, we see that the model exhibits two different administration approaches: either one with frequent infusions at lower dosage, or another one consisting in two repeated standard protocols (one infusion every 3 weeks for 4 cycles) with maximal doses. Both the values of t f t 1 I p (t)dt (tumor density) and Tox(t f , I p (t)) (toxicity) corresponding to the choice ∆ t = 3 weeks are higher than the ones obtained by imposing ∆ t = 1 week (Table 6). Similar results hold for the reference value of k AC (Tables 4 and 5), although the difference is lower. However, we cannot conclude that the optimal protocol with ∆ t = 3 weeks is less suitable than the one with ∆ t = 1 week, since our results need to be validated from the clinical point of view and our analysis does not consider side effects and practical constraints. On the one hand, we do not have any information about close-range ipilimumab infusions, and therefore the numerical results obtained with ∆ t = 1 week could not be reliable, since a patient could be affected by an increased toxicity depending on frequent infusions, which are not taken into account. Furthermore, we neglect the cost of the drug administration, which is, however, a crucial information to determine the feasibility of a therapy. For instance, comparing the mono-therapy protocols (rows 1 and 3 in Table 6), the one associated to ∆ t = 1 week consists in administering 12 infusions of ipilimumab, while the other one with ∆ t = 3 weeks consists in 8 infusions. As every infusion of ipilimumab costs around $30,000 [32], the predicted optimal administration protocol with ∆ t = 1 week, in reality, could not be feasible.
Our results point out the synergy between ipilimumab and sipuleucel-T (see Appendix B). By looking at the optimal protocols, the vaccine allows one to reduce the anti-CTLA4 dosage while reaching the goal of tumor eradication. The vaccine administered alone shows poor performances, while, coupled with the ipilimumab, its effectiveness considerably improves. The limited effect of the vaccine has been observed in other clinical studies [31], which highlight a few months increase in life expectancy, but not a PSA reduction. This could suggest to administer the dendritic cell vaccine in combination with other immunotherapies. Moreover, the results showed in Tables A1 and A3 predict a stronger synergy for lower values of k AC . This seems indicate that a combination therapy could provide best performances for those patients who do not have a good response to ipilimumab mono-therapy. The potential synergy between ipilimumab and sipuleucel-T has also been observed in clinical studies, such as [33], which report an increase of the median survival in patients with metastaticprogressive prostate cancer. A phase 1 clinical trial is currently investigating the efficacy of this combination immunotherapy in patients with advanced prostate cancer [34]. Thought a sensitivity analysis on a toxicity function parameter, we also investigate the reliability of the model prediction, by evaluating if the model outcomes can be highly influenced by the uncertainty in computing drug toxicity. Our analysis indicates that the computed optimal administration protocols seems not be dramatically affected by parameter perturbation, by highlighting a general robustness of the analysis herein reported.
In conclusion, the presented work is a theoretical analysis of the effect of ipilimumab under optimal protocols of administration. Additional efforts are needed to make this approach suitable to clinical applications. The definition we used for the toxicity function could be more accurate if more detailed experimental data were available. For example, we assumed that the resolution times of AEs of ipilimumab in PCa can be compared with the ones registered in melanoma. However, the median age of melanoma patients is usually lower than the one of PCa patients, and this could affect the resolution times. Moreover, the resolution time corresponding to the single infusion of ipilimumab was supposed to be the same as for the intermittent administration protocol, due to the lack of this information in the study by Small et al. [22]. When patient's treatment time series and associated clinical data will be publicly available, it would be also interesting to generate an in silico population of patients to compare the outcomes of the optimal control problem with clinical trial results, in terms of individual patient responses. It is important to note that our results depend on the initial state used for the model simulations. In the perspective of proposing a personalized protocol, the initial state should be set in accordance with the patient condition. The model can also be calibrated to patient-specific parameters, such as the tumor proliferation rate. In this perspective, the model outcomes can be used not only to find the optimal protocol for a personalized therapy, but also to evaluate if the patient responds well to the treatment. Indeed, by comparing the model predicted tumor progression with the patient prostate cancer evolution, the model could indicate, in almost real time, whether the treatment is working as expected, and therefore could be helpful in supporting clinical decisions.

Conclusions
In this paper, we employ a previously developed mathematical model of castrationresistant prostate cancer, in order to determine optimal protocols for the administration of the drug ipilimumab. This work provides a theoretical study about the efficacy of ipilimumab, working alone or in combination with the dendritic vaccine sipuleucel-T. Our model predicts that ipilmumab is potentially able to eradicate prostate cancer. Moreover, our results highlight a synergy between ipilimumab and sipuleucel-T, which allows for a reduction of the amount of ipilimumab administered, with a stronger effect on tumor reduction. Other efforts are necessary to employ our methods in a clinical environment, but we think that this theoretical work can help in understanding the potential effectiveness of the drug ipilimumab. Moreover, we hope that our results could encourage the scientific community to further investigate the role of immunotherapy agents in the treatment of advanced prostate cancer.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Toxicity Function Parameter Estimation
In order to determine the values of the parameters a, b and n of the Function (4), we implement a minimum problem aiming at reducing the squared distance between the toxicity coefficients and the function outputs: where Tox ] is a vector containing the ipilimumab toxicity coefficients in Table 2, and Tox([a, b, n], I p ) is the function defined in the Equation (4) By imposing ∂G ([a,b,n],I p ) ∂a ≥ 0, we obtain the optimal value: (A3) Therefore, the new minimum problem becomes: The parameter b and n have been estimated, respectively, between real and natural numbers in [1, 10 30 ] and [1,10], and the corresponding values are listed in Table 3.

Appendix B. Synergy
To investigate the synergy between the immunotherapies included in the model, we consider the Bliss combination index [35,36], which has been successfully employed to analyze the drug synergies in other studies [13,14]. This index compares the efficacy of the combination therapy, with the one obtained by considering the two drugs as not interacting, computed by probability laws. To do this, we define the efficacy of a treatment j as the resulting tumor reduction due to a therapy in comparison with the untreated case (T sx ): where j is the therapy, while t f = 5 years, is the time at which we evaluate the tumor reduction. The value of E(j) is in [0,1], where 0 means that the treatment j does not lead to a tumor reduction, while 1 indicates that T j (t f ) = 0, i.e., the cancer is eradicated. Therefore, we define a slightly modified Bliss combination index as: where E(AC, V) is the efficacy of the combination therapy including anti-CTLA4 and vaccine, while E(AC, 0) and E(0, V) are the efficacies of the two mono-therapies. If σ(AC, V) > 0, this means that the two drugs are synergistic, otherwise, if σ(AC, V) < 0, the two drugs are antagonistic.
Since the synergy defined in Equation (A6) depends on the treatment efficacies, we evaluate the combination indexes by varying the ipilimumab dose (Table A1). The synergy index is always positive, therefore the two drugs are predicted to be synergistic, but it first increases with ipilimumab dose, by reaching its maximal value for infusions of 3 mg/kg, and then starts decreasing.
This behavior could be easily explained by looking at the treatment efficacies, defined by Equation (A5). Indeed, since the model predicts a low effectiveness of the vaccine administered as mono-therapy, with E(0, V) = 1.98 · 10 −9 , the denominator of the Equation (A6) E(AC, 0) + E(0, V) − E(AC, 0)E(0, V) ≈ E(AC, 0). Therefore, the synergy index becomes a comparison between the the efficacy of the combination-therapy and the ipilimumab mono-therapy. Table A2 summarizes the model predicted efficacies of the ipilimumab mono-therapy according to the different doses. We observe that increasing the dose the estimated treatment efficacy tends to 1, which means that PCa after 5 years is close to zero, and therefore the addition of the dendritic cell vaccine does not lead to a significant improvement. Table A1. Index of synergy between sipuleucel-T (sip-T in the table) and ipilimumab for different ipilimumab doses, computed by the Formula (A6). The two drugs are administered by following the standard procedures. Anti-CTLA4 is administered every 3 weeks, for 4 cycles, while the vaccine is performed by infusing three doses of 50 million of cell every 2 weeks.  These results are influenced by the estimate of k AC , therefore we further analyze the synergy by considering a reduced value of this parameter. Table A3 shows the synergy indexes between sipuleucel-T and ipilimumab and the model predicted efficacies for the anti-CTLA4 mono-therapy for different ipilimumab doses. Table A3. Index of synergy between sipuleucel-T (sip-T in the table) and ipilimumab and anti-CTLA4 mono-therapy efficacy for different ipilimumab doses and an halved value of k AC . The coefficients of synergy and efficacy are computed by the Equations (A5) and (A6), respectively. The two drugs are administered by following the standard procedures. Anti-CTLA4 is administered every 3 weks, for 4 cycles, while the vaccine is performed by infusing three doses of 50 million of cell every 2 weeks. By comparing Tables A1 and A3, we observe that the model predicts a stronger synergy with a halved value of k AC than with the estimated parameter.  Figure A1. Comparison between administration protocols obtained with different values of the parameter k AC . The dynamics refer to the ipilimumab mono-therapy as result of the optimization problem OP2, with ∆ t = 3 weeks. Blue lines represent the model outcomes with the estimate value of k AC (the same presented in figure 6, blue line), while orange lines represent the one with k AC 2 . The ipilimumab threshold has been highlighted with a black dotted line. Figure A1. Comparison between administration protocols obtained with different values of the parameter k AC . The dynamics refer to the ipilimumab mono-therapy as result of the optimization problem OP2, with ∆ t = 3 weeks. Blue lines represent the model outcomes with the estimate value of k AC (the same presented in Figure 6, blue line), while orange lines represent the one with k AC 2 . The ipilimumab threshold has been highlighted with a black dotted line.  Figure A2. Comparison between administration protocols obtained with different values of the parameter k AC . The dynamics refer to the ipilimumab and sipuleucel-T combination therapy as result of the optimization problem OP2, with ∆ t = 1 week. Blue lines represent the model outcomes with the estimate value of k AC (the same presented in figure 5, orange line), while orange lines represent the one with k AC 2 . The ipilimumab threshold has been highlighted with a black dotted line. k Figure A2. Comparison between administration protocols obtained with different values of the parameter k AC . The dynamics refer to the ipilimumab and sipuleucel-T combination therapy as result of the optimization problem OP2, with ∆ t = 1 week. Blue lines represent the model outcomes with the estimate value of k AC (the same presented in Figure 5, orange line), while orange lines represent the one with k AC 2 . The ipilimumab threshold has been highlighted with a black dotted line.
parameter k AC . The dynamics refer to the ipilimumab and sipuleucel-T combination therapy as result of the optimization problem OP2, with ∆ t = 1 week. Blue lines represent the model outcomes with the estimate value of k AC (the same presented in figure 5, orange line), while orange lines represent the one with k AC 2 . The ipilimumab threshold has been highlighted with a black dotted line.  Figure A3. Comparison between administration protocols obtained with different values of the parameter k AC . The dynamics refer to the ipilimumab and sipuleucel-T combination therapy as result of the optimization problem OP2, with ∆ t = 3 weeks. Blue lines represent the model outcomes with the estimate value of k AC (the same presented in figure 6, orange line), while orange lines represent the one with k AC 2 . The ipilimumab threshold has been highlighted with a black dotted line. Figure A3. Comparison between administration protocols obtained with different values of the parameter k AC . The dynamics refer to the ipilimumab and sipuleucel-T combination therapy as result of the optimization problem OP2, with ∆ t = 3 weeks. Blue lines represent the model outcomes with the estimate value of k AC (the same presented in Figure 6, orange line), while orange lines represent the one with k AC 2 . The ipilimumab threshold has been highlighted with a black dotted line.   [15] ρ D Source of dendritic cells 0.00686 billion of cells day [15] µ D Death rate of dendritic cells 0.14 day −1 [37] λ Ipilimumab death rate 0.055 day −1 [15] i 0 Baseline level of IL-2 0.00299 ng mL [15] e Maximal activation rate of IL-2 by CTL 5000 ng mL·day·billion of cells [37] s

Appendix D. Parameters Estimates
Tumor cells saturation level for CTL stimulation 10 billion of cells [37]